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Abstract 

The method of self-similar factor approximants is shown to be very convenient for 
solving different evolution equations and boundary-value problems typical of physical 
applications. The method is general and simple, being a straightforward two-step 
procedure. First, the solution to an equation is represented as an asymptotic series 
in powers of a variable. Second, the series are summed by means of the self-similar 
factor approximants. The obtained expressions provide highly accurate approximate 
solutions to the considered equations. In some cases, it is even possible to reconstruct 
exact solutions for the whole region of variables, starting from asymptotic series for 
small variables. This can become possible even when the solution is a transcendental 
function. The method is shown to be more simple and accurate than different variants 
of perturbation theory with respect to small parameters, being applicable even when 
these parameters are large. The generality and accuracy of the method are illustrated 
by a number of evolution equations as well as boundary value problems. 
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1 Introduction 



Differential equations appear in numerous problems of physics, applied mathematics [1-3] 
and many other branches of natural as well as social sciences (see, e.g., [4-7]). In the majority 
of cases, these equations are nonlinear and cannot be solved exactly, allowing for exact 
solutions only in a few exceptional instances. Then, in order to obtain approximate analytical 
solutions, one resorts to perturbation theory in powers of some parameters assumed to be 
small [1-3]. The resulting expressions are usually rather cumbersome and are difficult to 
analyze. Also, they form asymptotic series that are useful only for very small expansion 
parameters. 

The validity of perturbative series can be extended to the finite values of parameters by 
reorganizing them with the help of the optimized perturbation theory [8]. The basic idea 
of this theory is to include in the initial approximation a set of auxiliary parameters which 
are transformed, at each step of perturbation theory, into control functions governing the 
series convergence [8]. The optimized perturbation theory has been successfully applied to 
a great variety of problems, providing rather accurate approximations (see review article [9] 
and references therein). It has also been applied to solving differential equations [10-17]. 
However, the weak point of this approach to finding the solutions of differential equations is 
that the optimization procedure results in very complicated equations for control functions, 
which are to be solved numerically. It is practically always much easier to solve the given 
differential equations numerically than to deal with the cumbersome optimization equations 
for control functions, anyway requiring numerical solution. This is why this approach, though 
being very useful for many other problems [9], has not found wide practical use for solving 
differential equations. 

Another method of constructing approximate solutions is based on the self-similar ap- 
proximation theory [18-26]. Then the solutions to differential equations can be represented 
in the form of self-similar exponential approximants or self-similar root approximants [27- 
34]. These approximants represent well those functions whose behavior at large variables is 
known to be either exponential or power-law, respectively. 

In the present paper, we advocate a novel approach to constructing approximate solu- 
tions of differential equations. This approach is based on the use of self-similar factor ap- 
proximants [35-37]. The mathematical derivation of the latter also rests on the self-similar 
approximation theory [18-26], but the structure of these approximants is rather different 
from the exponential and root approximants [27-34]. The structure of the self-similar fac- 
tor approximants reminds that of thermodynamic characteristics near critical points [38,39]. 
This is why it has been natural to apply, first, these approximants to the description of criti- 
cal phenomena [35-37,40]. It was shown that these approximants allow us a straightforward 
and simple determination of critical points and critical indices, agreeing well with the results 
of the most complicated numerical techniques, whose description can be found in articles 
[41-46] and books [47-49]. 

The self-similar factor approximants make it possible to define an effective sum of diver- 
gent series. Initially, these approximants were introduced [35-37] for summing the partial 
series of even orders, while the summation of odd-order series was not defined. Recently, the 
method was completed by defining the factor approximants of odd orders [40,50]. Now we 
have in hands a general and uniquely prescribed procedure for constructing the self-similar 
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factor approximants of arbitrary orders. We show below that this procedure can be employed 
for finding very accurate approximate solutions to differential equations. In some cases, the 
method gives exact solutions, if these exist. The approach is very general, being applicable 
to linear as well as to nonlinear equations, and to initial-value as well as to the boundary- 
value problems. The principal difference of the method from other perturbation theories 
[1-3] is that we, first, represent the sought solutions as an asymptotic series in powers of the 
equation variable, but not in powers of a parameter. And, second, we extrapolate the given 
series to the whole range of the variable by means of the self-similar factor approximants. 
The advantage of the method is its extreme simplicity combined with the high accuracy of 
the obtained solutions. 

2 Self-similar factor approximants 

Here we describe the general procedure of constructing self-similar factor approximants as 
solutions to differential equations. We keep in mind ordinary differential equations, though 
the procedure can be generalized to partial differential equations. It would be unreasonable 
to plunge from the very beginning to complicated matters, but the principal idea should, 
first, be illustrated by not too complex problems. In line with this, we consider the case, 
where the solution is a real function of a real variable. Thus, we assume that we are interested 
in finding a solution to a differential equation, which is a real function f(x) of a real variable 
x. First, we are looking for the solution in the region of asymptotically small x — > 0, where 
it can be represented as an asymptotic series 

f{x)~f (x) ]Ta n x" (x-0), (1) 

n 

where fo(x) is a given function. It is worth emphasizing that the expansion variable x is not 
compulsory the given independent variable, but can be a function of it, which means that 
one can always accomplish the change of variables and look for the expansion in powers of 
the new variable. For instance, if the initial variable entering an equation is t, we are not 
obliged to look for an expansion necessarily in powers of t, but we can make a change of the 
variable, introducing x — x(t) and to study an expansion in powers of x(t). In particular, 
x(t) could be t a . If a is not an integer, we get the Puiseux expansion. In general, one 
can consider asymptotic expansions of the form J2n a nX n , with functions x n from the given 
comparison scale [51]. 

Suppose, series (1) is limited by a finite sum of fc-th order 

f k (x) = f (x) £ <>nX n , (2) 
n=0 

with k — 0,1,2, . . .. Without the loss of generality, we may set a — 1, since if a would not 
be 1, we could include it in fo(x). 

The self-similar factor approximants of even orders k = 2p, with p — 1,2, . . ., are defined 
[35-37] as 

fi P (x)=f (x)f[(l + A i x) nt . (3) 

i=l 



3 



The parameters A-i and rii are obtained from expanding Eq. (3) in powers of x and com- 
paring the resulting expansion with sum (2) by the accuracy-through-order procedure. The 
corresponding equations for aj and rij can be written as 



Y j n l A n l =B ri 



[n 



1,2, 



i=l 



where 



B r , 



- — lim - — In a m x r 

(n-1)! ~o dx" ^ 



(4) 



(5) 



Equations (4), whose number is 2p, define all 2p unknowns Ai and rij. The factor approxi- 
mants (3) extrapolate expansion (2), valid for small x — > 0, to the region of finite x. 

For odd orders fc = 2p + 1, with p = 0, 1,2, . . ., the factor approximants are defined 
[40,50] in the form 



(6) 



i=i 



in which Ai = 1 and all other 2p + 1 unknown parameters are given by 2p + 1 equations 



p+i 

EM? 

1=1 



[n 



l,2,...,2p+l), 



(7) 



with _B n from Eq. (5). 

The index k, labelling the factor approximants fk(x), is determined by the order k of 
fk(x) in sum (2), which is used for calculating the parameters Ai and rij. It may happen that 
there exist additional conditions imposed on the behavior of the sought function f(x). Then 
the structure of f£(x) should include additional factors, whose parameters are chosen so 
that to satisfy the imposed conditions. In the case of differential equations, such additional 
restrictions are given by initial and boundary conditions. 

For example, the Cauchy initial value problems are supplemented by the set of initial 
conditions for the sought function and its derivatives 



/(0) = fo , 



d m f(x) 
dx m 



(m) 



(8) 



x=0 



where m — 1,2, Therefore, the factor approximants are constructed so that to satisfy 

the same initial conditions 



tf(0) = /o, 



d m r k {x) 

dx m 



Jo 



(m) 



(9) 



x=0 



for all orders of k. 

In the case of a boundary value problem, for a function f(x) on an interval of x G [x±, x 2 ], 
one has the boundary conditions 



f(xi) = fi , f(x 2 ) = h ■ 



(10) 
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Respectively, the factor approximants of each order have to satisfy the same boundary con- 
ditions 

KM = h , fi(x 2 ) = h ■ (11) 

When a boundary value problem is formulated for an infinite interval, the related bound- 
ary conditions are given in the form of the asymptotic behavior of the function f(x). For 
instance, let the asymptotic form of f(x) at large x be prescribed as 

f(x) ~ Bx? (x -> oo) , (12) 

where B ^ and (5 is any real quantity, including zero. Hence, the factor approximants 
fk(x) of any order must possess the same asymptotic form 

f* k {x)~BxP {x^oo). (13) 

Assume that the zero-order factor in Eqs. (1) and (2) behaves at large x as 

f (x) ~ cx a (x -> oo) , (14) 

with c 7^ and any real a. Then, in order that the factor approximants (3) and (6) would 
satisfy the asymptotic boundary condition (13), it should be that 

c]\A?=B, a + £> = /3. (15) 

i i 

Depending on whether initial-value or boundary-value problems are considered, we shall 
always require that the corresponding conditions be valid exactly. This distinguishes our 
approach from perturbation theories with respect to some parameters, when boundary con- 
ditions are usually satisfied only approximately [1-3]. 



3 Initial value problems 

Considering the Cauchy initial value problems, we shall concentrate our attention on the 
so-called singular problems, where a parameter enters in front of the highest derivative [1- 
3]. When one tries to apply perturbation theory with respect to the given parameter, one 
confronts nonuniformly valid expansions [1-3]. However in our approach, we start with an 
expansion not in powers of this parameter, but in powers of the equation variable. Then, 
using the extrapolation by self-similar factor approximants, we obtain accurate solutions for 
any value of the parameter. Moreover, in some cases, the factor approximants restore exact 
solutions. This fact is so unusual that we illustrate it in detail below. 

A. Linear singular problem 

Let us consider the singular initial value problem 
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with the initial conditions 



defining y — y(t) for t > 0. And let us pretend that we do not know its exact solution. 

Following our approach, we, first, look for the behavior of the solution y(t) at asymptot- 
ically small t — > 0, which is 

y(t) ~ 1 - + ^ + . . . (f -> 0) . (18) 

For the factor approximant of second order, we get 

y 2 *(t) = hm(l + A0- 1/A£ = e-*/ £ , (19) 

which is the exact solution of Eq. (16). Similarly, for the factor approximants of higher 
orders after k — 3,4,5 and so on, we find 

y* k (t) = (k > 3) , (20) 

that is, the exact solution of Eq. (16) under the initial conditions (17) and any real e. 

It is easy to notice that function (20) cannot in principle be expanded in terms of e. The 
functions of this type are called instanton solutions [47,52]. Note also that the exponential 
(20) is a transcendental function. 



B. Nonlinear singular problem 

The model problem of carrier transfer [2,3] is given by the nonlinear equation 

(ey + t)^ + y-l = (21) 

for y = y(t) and t > 0, with the initial condition 

2/(0) = 2 . (22) 

We again pretend that we are not aware of the solution to Eq. (21). If we try to invoke 
perturbation theory with respect to the parameter e, we come to a rather nontrivial singular 
problem. The method of strained coordinates or the method of asymptotic matching could 
be used, but with a necessity of dealing with very lengthy and cumbersome calculations [2,3]. 
We show here that this problem is easily solvable by our method. 

As is explained in Sec. 2, it is always possible to make a change of variables, which seems 
to be convenient. If we introduce the function z(x), such that 



then Eq. (21) simplifies to 



t . . 

z — x + y , x = - , (23) 

£ 



dz . , 

z- x - 1 = 0, 24 

ax 
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with the initial condition 

z(0) = 2 . (25) 
The asymptotic behavior of z(x) at small x is 

1 3 

z(x) ~2 + - x+ — x 2 + ... (x -> 0) . (26) 

Constructing the factor approximants, we find that, starting from the fourth order, the 
approximants z%(x) provide the exact solutions to Eq. (24), with the initial condition (25). 
Thus, in fourth order, we have 

zl(x) = 2(1 + A 1 x) ni (l + A 2 x) na . (27) 

From the accuracy-through-order procedure, with the given initial condition (35), we get 

A 1 = 1 (l - iV3) , A 2 = A{, m = n 2 = X - . (28) 

Thence Eq. (27) reduces to 

z\{x) = V4 + 2x + x 2 , (29) 

which is the exact solution of Eq. (24). Acting in the same way for higher-order approxi- 
mants, and using relations (23), we come to the factor approximants 



2t t 2 t 

y^ = f + j + -2 - - £ (*^ 4 ). (3°) 

giving the exact solution of Eq. (21) for orders k > 4, any real e, and all t > 0. 
C. Singular logistic equation 

The logistic equation and its variants are widely employed in the studies of population 
dynamics [53]. For the population p(t), being a function of time t, with the growth rate 1/e, 
the logistic equation writes as 

ef t =p(l-p). (31) 
It is given for t > 0, with the initial condition 

p(0) = Po . (32) 

If again we pretend not to know the solution of Eq. (31) and try to resort to perturbation 
theory in powers of the large growth rate, that is, small e, we confront a singular value 
problem. 

Wishing to apply our approach, it is reasonable, as has been discussed above, to choose 
the most convenient variables. The hint on what variables would be the most convenient 
comes from analyzing the dynamics in the vicinity of the stationary points. Equation (31) 
possesses one stable fixed point p(oo) = 1. The motion near this point is given by the 
asymptotic law 

p(t) ~ 1 + ce~ t/£ (t -> oo) . (33) 
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The latter suggests to choose as the variable 

x = e- t/£ . (34) 
Considering the population p(t) as a function of this variable (34), for the function 

y = p(t(x)) = y(x) , (35) 

we have the equation 

^ = 1/(^-1)^, (36) 

with the initial condition 

y(l)=Po- (37) 

The small x-expansion of y(x) is 

y(x)~Y,a n x n (x-0), (38) 

n 

for which Eq. (36) yields a n = a™. The lowest-order factor approximant, for which the initial 
condition can be satisfied, is the third-order approximant resulting in the expression 

yZ(x) = (l + A lX ) n \ (39) 

where 

Ax = l , m = -l . 

Po 

The same expression (39) follows for other approximants of higher orders k > 3. Returning 
to the population function p(t) by means of relations (34) and (35), we find 

aw = Be-* {k> - 3) - m 

This is the exact solution of the logistic equation (31), with the initial condition (32). 



4 Boundary value problems 

In the previous section, we have shown that the method of self-similar factor approximants 
allows us to reconstruct exact solutions of some differential equations representing initial 
value problems. This highly nontrivial fact occurs as well for some nonlinear boundary value 
problems. 



A. Kink soliton equation 

Nonlinear equations possessing soliton solutions are met in various problems of physics 
and applied mathematics. Let us consider, for example, the nonlinear Schrodinger equation 
describing the so-called y9 4 -model with particle mass \ je. The equation reads as 

e d 2 (p i 3 



2 dx 2 



+ <p-<p 6 = 0. (41) 
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Assume that the function ip — ip(x) satisfies the boundary conditions 

(p(-oo) = -1 , (p(oo) = 1 . (42) 

From Eq. (41) and the boundary conditions (42) it follows that ip(x) is an antisymmetric 
function, such that 

<p(-x) = -<p(x) , (43) 

and, therefore, 

<p(0) = . (44) 

To choose a convenient expansion variable, we again study the form of the function tp(x) 
in the vicinity of the stable stationary points </?(=Foo) = =Fl, where we have 

ip(x) ~ — 1 + aexp f —= x ) (x — > — oo) , 
( 2 \ 

(p(x) ~ 1 — aexp ■= x (x — > oo) , (45) 

with a being a real parameter. This immediately suggests to choose as a variable 

(2 \ 

z = exp —= x . (46) 
Vv^ J 

In terms of the latter variable, Eqs. (45) become 

<f(x) ~ — 1 + (z — > 0, rr — > — oo) , 

ip(x) ~ 1 — - (z — > oo, x — > 0) . (47) 
It is also convenient to introduce the positively defined function 

y = 2 + <p = y{z) , (48) 

satisfying the boundary conditions 

y(0) = 1 , y(l) = 2 . (49) 

Then Eq. (41) transforms into 

2z 2 ^ + 2z^- + Q-lly + Qy 2 -y 3 = . (50) 

dz z dz 

Looking for the form of y(z) at small z — > 0, we substitute the asymptotic expansion 

y{z)~Y. a nZ n ("0) (51) 

n 

into Eq. (50), from where it follows that 

ao = l, a n = 1 a\ (n>l). (52) 
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Constructing the factor approximants from expansion (51), we obtain the same expression 

vttz) = (* > 4) (53) 

2 + a±z 

for orders k > 4, with the constant a x = 2 found from the boundary conditions (49). Using 
relations (46) and (48), according to which 

<p(x)=y(z(x))-2, (54) 

we come to the expression 

^(s) = tanh ^ (fc>4). (55) 

This is the exact kink solution of the soliton Eq. (41). 

Recall that, if we would try to solve Eq. (41) by means of perturbation theory with 
respect to large mass, that is, small e, we would have to deal with a rather unpleasant and 
cumbersome boundary- layer problem [1-3]. 

B. Bell soliton equation 

The nonlinear Schrodinger equation for a negative mass — \je is 

5 ^ - V + ^ = 0. (56) 
The function if — if (x) satisfies the boundary conditions 

if(-oo) = , if(oo) = . (57) 
From these it follows that if{x) is a symmetric function, 

if(-x) = if(x) , (58) 

such that 

lim ^ = . (59) 
z->±o dx 

In the vicinity of the stationary points y?(=Foo) = we have 

if(x) ~ a exp | ij— x\ (x — > — oo) , 

if{x) ~ a exp | — ]J - x ] (x — > oo) . (60) 
Hence the appropriate convenient variable here is 

z = exp I J- x | . (CI ) 
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In terms of the latter, Eq. (60) simplify to 

(p(x) ~ az (z — > 0, x — > — oo) , 

</?(:r) ~ - (-2 — > oo, x — > oo) . (62) 

z 

For the function 

y = ¥>(z(*)) = y(^) (63) 

of variable (61), from Eq. (56), one gets the equation 

*_ + *_-,, + „ =0, (64) 

with the boundary conditions 

y(0) = , y(oo) = . (65) 
And property (59) takes the form 

lim ^ = . (66) 
In the asymptotic region of small z, the solution to Eq. (64) is 

y(z)^^a n z n (*-0), (67) 

n 

in which 

a 2n = (n = 0,1,2,...) , 

a 2n+ i= ("£f (68) 
The corresponding factor approximants take the same form 

«*M = rn^F (fc£3) (69) 

for orders > 3, with the parameter a± = 2\^2 given by the boundary condition (66). 
Resorting to relations (61) and (63), we obtain 

<p* k (x) = V2 sech (k > 3) , (70) 

which is the ea;aci solution of Eq. (56) describing a bell soliton. 
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5 Evaluation of approximation accuracy 



In the previous sections, we have considered several examples of initial-value and boundary- 
value problems, for which the self-similar factor approximants result in exact solutions. 
This fact, as such, that starting with approximations, one can get exact solutions, is highly 
nontrivial. This is why we have focused our attention on its thorough illustration. But, 
certainly, the more general situation is when an exact solution cannot be reconstructed, or 
just does not exist at all. How then could we evaluate the accuracy of our approximate 
solutions? 



A. Solutions defects and errors 

Suppose a differential equation, that can be represented in the operator form as 

E[y(x)]=0, 

defines a function y(x) of a variable x in the interval [rci,^]. The latter can be finite or 
infinite. Let us find an approximate solution yl(x) of the A;-th order. In the theory of 
differential equations [54,55], one characterizes the accuracy of approximate solutions in two 
ways, by calculating the solution defects and solution errors. 
The solution defect of y k (x) is 

D[y* k (x)] = | E[y* k (x)] | . (71) 

This is a local characteristic of an approximate solution yl(x), showing to what extent the 
considered solution y k (x) does not satisfy the given equation E[y(x)] = 0. For the exact 
solution y(x), the solution defect is, evidently, zero. 

Varying x in the whole interval [xi,^], one defines the maximal solution defect 

D[y* k ] = sup D[yt(x)\ , (72) 

X 

which is the global characteristic of the validity of the approximate solution yl(x), with 
respect to the given equation E[y(x)] = 0, in the whole interval of the variable x G [xi,a; 2 ]. 
The solution error of y k (x) is defined as 

A[y* k (x)} = \y* k (x) - y(x)\ . (73) 

This local characteristic shows how much the approximate solution yl(x) deviates from the 
exact solution y(x) at each x. 
The maximal solution error 

A[yl\ = sup A[yl(x)\ (74) 

X 

gives the global characteristic of the maximal deviation of y k (x) from y(x) in the total interval 
[x 1 ,x 2 \. 

Both the solution defects and solution errors characterize the accuracy of approximate 
solutions. The explicit relation between these characteristics depends on the explicit form 
of the given differential equation. Also [54,55], one can define the error-to-defect ratio 
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B. Boundary layer problem 



To illustrate the definitions of the previous subsection in the case of the factor approxi- 
mants, let us consider the boundary-layer problem given by the equation 

, d y n c?fi\ 

£ d^ + x T x - xy = ^ (76) 

where the function y = y(x) is defined for x in the interval < x < 1, with the boundary 
conditions 

y(0) = , y{l) = e . (77) 

Equation (76) is an example of a nontrivial boundary-layer problem [2,3] for which the 
conventional matching technique contains not only powers of e, but also powers of yfe, and 
Ine, which makes the matching rather complicated. Also, perturbation theories in terms of e 
do not allow for the exact validity of the boundary conditions (77), but the latter are satisfied 
only approximately, for small e C 1, making such theories [2,3,12] inapplicable for large e. 
Contrary to this, the method of factor approximants, being very simple, gives approximate 
solutions for arbitrary e, with high accuracy, and satisfies the boundary conditions exactly 
for any e. 

It is convenient, first, to redefine the sought function as 

z = ye-* = z(x) , (78) 

in order that the latter be varying between zero and one. For function (78), problem (76) 
transforms to 

d 2 z ( x\ dz ,_„, 

5? + ( 2 + 7)^ = ' (79) 



with the boundary conditions 



z(0) = , z(l) = 1 . (80) 



Then we follow the standard procedure, by deriving an expansion of z(x) for asymptotically 
small x — > 0, and by constructing the factor approximants based on the latter expansion. 
Finally, using relation (78), we obtain the solutions of Eq. (76) as factor approximants yl(x). 

The solution defects (71) for several factor approximants yl(x) are shown in Fig. 1, which 
demonstrates good uniform convergence of the method. Figure 2 shows the solution errors 
(73) for the same factor approximants. To better investigate the relation between different 
accuracy characteristics, we present the maximal solution defects (72), maximal solution 
error (74), and the error-to-defect ratio (75) in Tables 1 and 2. For small e, the numerical 
convergence is slower than for large e. Therefore, in Table 1 for e — 0.1, we present the 
accuracy of higher approximants up to k = 17. Numerical convergence for e > 1 is so fast 
that it is sufficient to consider the approximants up to k = 7, as in Table 2. As is seen, 
the maximal solution errors are much smaller than the maximal solution defects. Hence the 
latter can serve as an upper bound for the former. 



13 



6 Nonlinear radial equations 

Now we shall illustrate the power of the method for some physically motivated nonlinear 
equations in radial variables. 

A. Gross-Pitaevskii equation 

This is the nonlinear Schrodinger equation applied to diferent Bose-condensed systems 
(see Refs. [56-62]). We consider here the variant of this equation describing vortices in Bose 
systems. Analogous equations describe also vortices in superfluids, superconductors, and 
Higgs fields. The equation, in dimensionless units, reads as 

d 2 ip ^ 1 d(p ^ + ip (/9 3 — (81) 
dr 2 r dr r 2 

defining a function if — ip(r) of the radial variable r > 0. The boundary conditions are 

<p(0) = , y?(oo) = 1 . (82) 
The asymptotic behavior of the solution to Eq. (81) at small r is 

V?(r) ^ 5> 2n+1 r 2 ™ +1 (r-0), (83) 

n 

where 

1 1 + 8a 2 

«3 = - g ai , a 5 = a, , 

1 + 80a? 1 + 656a? + 1152a? 

a 7 = ai , ao = ai , 

7 9216 737280 

and so on, with a\ to be defined later from the second of the boundary conditions (82). For 
convenience, expansion (83) can be rewritten as 

V9(r)~ ai r $> n r 2 " (r 2 -> o) , (84) 

n 

where b n = a 2n+ i/a 1 . Expression (84) shows that actually, we have an expansion in powers 
of r 2 . The corresponding factor approximants have the form 

^w=^n(i+v) ni . (85) 

i 

The boundary condition 

<pl(0) = (86) 

for Eq. (85) does not impose additional constraints, being always valid. And the boundary 
condition 

rt(oo) = 1 (87) 

imposes on Eq. (85) two constraints 

c k l[A? = l, l + 2^n, = 0. (88) 
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In labelling the factor approximants, we associate the order k with the number of terms 
in expansion (84) with respect to r 2 . Then for the second-order approximant ip 2 , defined by 
Eq. (85), we have 

c 2 = 0.518840 , A 1 = 1 , A 2 = 0.287401 , 

n x = 0.026243 , n 2 = -0.526243 . 
For the third-order approximant ip 3 (r), we get 

c 3 = 0.585667 , A l = 0.107803 , A 2 = 0.139245 , 

ni = 1.761220 , n 2 = -2.261220 . 
The fourth-order approximant (p%(r) is defined by 

c 4 = 0.585331, A 1 = 1 , A 2 = 0.117323 , 

A 3 = 0.130445 , m = -0.000031 , n 2 = 4.553526 , n 3 = -5.053495 . 
For the fifth-order approximant <fl(r), we find 

c 5 = 0.583142, Ai = 0.158576 , A 2 = 0.045537 + i 0.011910 , A 3 = A* 2 , 

n x = -0.994835 , n 2 = 0.247422 - i 0.429221 , n 3 = n* 2 . 

Similarly, any higher-order approximant can be defined following the standard procedure 
of Sec. 2. Numerical convergence is quite fast, which is shown in Table 3. Note that 
the approximate solutions, for the same Eq. (81), found by employing the self-similar root 
approximants [27], are less accurate than the self-similar factor approximants, considered 
here. To compare the accuracy of the factor approximants tp* k {r) and the root approximants 
Rl{r), given in Appendix A, we present in Table 3 both the maximal solution defects -D 
as well as D[Rl\. As is seen, the factor approximants are two orders more accurate than the 
root approximants. 

B. Stokes-Oseen equation 

A simplified spherically symmetric variant of the Stokes-Oseen equation [1-3] can be 
written in the form 

d 2 u ^ 2 du ^ du q 
dr 2 r dr dr 

where u = u(r) and r > 1. This equation describes the viscous flow past a sphere of unit 
radius, with e playing the role of the Reynolds number. The boundary conditions are 

«(1) = , u(oo) = 1 . (90) 

When considering this equation by means of perturbation theory with respect to the Reynolds 
number e, one confronts a very delicate singular boundary-layer problem [2,3], with compli- 
cated expansions and matching, involving unexpected orders such as eln(l/e). But in our 
approach, the problem is easily solvable. 



15 



It is convenient, first, to change the variable to 



The function 
obeys the equation 

with the boundary conditions 



x = r-l. (91) 
V = u(r(x)) = y(x) (92) 



— + 2 < ^L + ey < ^L = Q (93) 
dx 2 1 + x dx dx 



y(0) = , y(oo) = 1 . (94) 



Then, as usual, we construct the factor approximants yl(x) and, using relations (91) and 
(92), return to the factor approximants u* k (r). The maximal defects of the latter are given 
in Table 4 for different e. For small Reynolds numbers e, the numerical convergence is very 
fast. For e ~ 1, convergence is a little slower, though again it is easy to reach high accuracy. 
It is feasible to reach quite good accuracy even for large e » 1, although this requires to 
construct the factor approximants of higher orders k ~ 10. 

C. Strongly-singular problem 

It is instructive to analyze the equation 

d 2 u 1 du / du\ 2 du 
dr 2 r dr \dr J dr 

defining the function u(r) for r > 1, with the same boundary conditions 

u(l) = , u(oo) = 1 , 

as in Eq. (90). If this equation is treated by perturbation theory [2,3] with respect to 
e, it becomes a terribly complicated problem. The matching procedure for that equation 
is notoriously difficult, because an infinite number of terms, with respect to e, must be 
calculated before even the leading order could be explicitly matched [2,3]. 

In the method of the self-similar factor approximants, the procedure is straightforward, 
as is described in Sec. 2. We follow the same way as in the previous subsection. A little 
difference from the previous problem (89) is that for Eq. (95) not all even-order approximants 
exist, but when they exist, they are close to the odd-order approximants, similarly to the 
case of Eq. (89). The odd-order approximants exist for all orders we have checked and 
their maximal defects are presented in Table 5 for different e, demonstrating good numerical 
convergence and the possibility of reaching quite good accuracy. 

7 Conclusion 

We have shown that the method of self-similar factor approximants [35-37,40,50] can serve 
as a powerful tool for constructing approximate solutions to ordinary differential equations. 
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These can be the initial-value and boundary-value problems, linear as well as nonlinear 
differential equations. 

A very nontrivial fact is that in some cases factor approximants automatically reconstruct 
exact solutions even to nonlinear equations, provided such solutions exist. This happens, 
when the latter solutions can be represented, with the appropriate change of variables, as 
factor approximants, that is, when the solutions pertain to the class of exactly reproduceable 
functions [35,50]. 

The procedure of constructing the factor approximants is very simple, which is described 
in Sec. 2. This procedure also is uniquely defined. If the asymptotic expansion of a solution 
has the form of Eq. (1), then all factor approximants enjoy the same structure 

r k (x) = f (x) U(i + Axp, 

i 

whose parameters are defined by the accuracy-through-order procedure. 

The quantities A; L and rii can be functions of all other parameters entering the considered 
equation. This is why the use of the factor approximants is not limited to the case of small 
parameters, but is valid for any values of the latter, yielding approximate solutions of high 
accuracy. 

As has been proved in our previous publications [35-37,40,50,63] the self-similar approx- 
imants guarantee an essentially higher accuracy than Pade approximants [64]. This fact 
can be easily understood, if one notices that Pade approximants have the structure which is 
just a particular variant of the factor approximants. In addition, the factor approximants, 
contrary to Pade approximants, enjoy a great advantage of being uniquely defined. 

Finally, since the quantities and rij in the expression for the factor approximants can 
be functions of other parameters, they also can be functions of other variables, entering the 
equations, and even functions of external fields. This could open the way for extending the 
suggested technique to solving partial and stochastic differential equations [65]. The latter, 
however, is the subject for future investigations. 
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Appendix A 



In Sec. 6, subsection 6.1, the factor-approximant solutions <fl(r) to the Gross-Pitaevskii 
equation (81) are compared with the root approximant solutions R* k (j')- In Table 3, the 
maximal defects of both types of approximants are given. 

The self-similar root approximants [26,27,29,30] are constructed from the asymptotic 
expansion 

/ x n 1 -2 9 _ 4 161 _ 6 
yy > ~ 2 8 16 

for large r — > oo, which is an expansion in powers of r~ 2 . The first several root approximants 
are 



Rl{r) = 



r 



1 „\- 1/2 



r 2 



2 V 4 

1 1 \ - 1 / 4 

J- 1 A \ 



V2 

The maximal solution defect D[Rl] happens at r ~ 1 and is shown in Table 3. Though the 
accuracy of the root approximants is not bad, it is two orders lower than that of the factor 
approximants. 
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Figure Captions 



Fig. 1. Solution defects D[y* k (x)} of several factor approximants y* k (x) with e — 1, for 
k = 4 (solid line), = 5 (dashed line), k = 6 (dotted line), and = 7 (dashed-dotted line), 
forEq. (76). 

Fig. 2. Solution errors A[y|(x)] of the same factor approximants, as in Fig. 1, for 
the same orders k = 4 (solid line), k = 5 (dashed line), k = Q (dotted line), and k = 7 
(dashed-dotted line), for Eq. (76). 



22 



Table Captions 

Table 1. Maximal solution defects D[yl], maximal solution errors A[yj£], and the error- 
to-defect ratio 5[y k ] for Eq. (76), with e = 0.1. 

Table 2. Maximal solution defects D[yl\, maximal solution errors A[yj£], and the error- 
to-defect ratio 5[yl\ for Eq. (76), with e — 1 and e = 10. 

Table 3. Maximal solution defects -D[y?£] of the factor approximants (p%(r) and maximal 
solution defects D[R* k ] of the root approximants R* k (r) for Eq. (81). 

Table 4. Maximal solution defects of the factor approximants u* k (r) for different 

Reynolds numbers e, in the case of Eq. (89). 

Table 5. Maximal solution defects D[u* k ] of the factor approximants u* k (r) for different 
e, in the case of the strongly singular Eq. (95). 
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Figure 1: Solution defects D[y%.(x)] of several factor approximants y%{x) with e — 1, for 
k = 4 (solid line), k = 5 (dashed line), A; = 6 (dotted line), and k = 7 (dashed-dotted line), 
forEq. (76). 
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Figure 2: Solution errors A[y^(x)] of the same factor approximants, as in Fig. 1, for the same 
orders k = 4 (solid line), k — 5 (dashed line), k = 6 (dotted line), and k = 7 (dashed-dotted 
line), for Eq. (76). 
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k 


D[yt\ 






8 


1.8 


0.14 


0.078 


9 


0.58 


0.051 


0.088 


10 


0.20 


0.018 


0.090 


11 


0.68 


0.038 


0.056 


12 


0.28 


0.015 


0.054 


13 


0.090 


0.0048 


0.053 


14 


0.024 


0.0014 


0.058 


15 


0.057 


0.0022 


0.039 


16 


0.032 


0.0012 


0.038 


17 


0.0083 


0.00028 


0.034 



Table 1: Maximal solution defects D[yl\, maximal solution errors A[y£], and the error-to- 
defect ratio S[yl] for Eq. (76), with e — 0.1. 







- 1 






: - 10 




k 


D[yt\ 


A|j/2] 




D[yl\ 


m\ 


S[yl\ 


4 


2.0 


0.062 


0.031 


1.58 


0.0051 


0.0032 


5 


0.42 


0.0094 


0.022 


0.059 


0.00013 


0.0022 


6 


0.20 


0.0040 


0.020 


0.025 


0.000046 


0.018 


7 


0.055 


0.00064 


0.012 


0.0032 


3.8-10- 6 


0.012 



Table 2: Maximal solution defects D[yl], maximal solution errors A[y£], and the error-to- 
defect ratio 5[yl] for Eq. (76), with e — 1 and e = 10. 



26 



k 


D[y* k \ 


D[RV\ 


2 


0.12 


0.14 


3 


0.017 


0.25 


4 


0.015 


0.10 


5 


0.0020 


0.11 


6 


0.0018 


0.10 



Table 3: Maximal solution defects -D of the factor approximants (f* k {r) and maximal 
solution defects D[R* k ] of the root approximants R k (r) for Eq. (81). 



k 


e = 0.1 


- 1 


: - 10 


4 


0.86 


0.094 


3.1 


5 


0.0035 


0.015 


12 


6 


0.00068 


0.024 


40 


7 


0.00015 


0.023 


12 


8 


0.00010 


0.00095 


0.32 


9 


0.000056 


0.00095 


0.15 


10 


0.000056 


0.0015 


0.15 


11 


0.000026 


0.0015 


0.23 


12 


0.000024 


0.0015 


0.011 



Table 4: Maximal solution defects D[u* h ] of the factor approximants u* k (r) for different 
Reynolds numbers e, in the case of Eq. (89) 



k 


e = 0.1 


: = 1 


= 10 


3 


0.0035 


0.062 


1.4 


5 


0.00056 


0.080 


0.44 


7 


0.00060 


0.011 


0.54 


9 


0.00035 


0.0017 


0.32 


11 


0.00015 


0.0017 


0.078 



Table 5: Maximal solution defects D[u* k ] of the factor approximants u* k (r) for different e, in 
the case of the strongly singular Eq. (95). 
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